
#pragma once

#include <stdexcept>
#include <string>
#include <stdio.h>
#include <math.h>
#include "kernel_defs.hpp"
#include "maths_functions.hpp"
/*//commented by Harsha start
// Note: coef should be e^-2*pi*i*delay
template<gpu_size_t BLOCK_SIZE, typename Complex, typename Real>
__global__
void phase_shift_kernel(const Complex* d_in,
                        Complex*       d_out,
                        Real           coef,
                        gpu_size_t     count)
{
	gpu_size_t thread_idx = threadIdx.x + blockIdx.x*BLOCK_SIZE;
	gpu_size_t thread_count = BLOCK_SIZE * gridDim.x;
	
    // Iterate over values in the input
	for( gpu_size_t idx=thread_idx; idx<count; idx+=thread_count ) {
		Complex a = d_in[idx];
		
		Complex phase;
		// TODO: Check if pre-computing these and loading from mem is faster
		__sincosf(coef*idx, &phase.y, &phase.x);
		
		Complex result;
		result.x = a.x*phase.x + a.y*phase.y;
		result.y = a.y*phase.x - a.x*phase.y;
		d_out[idx] = result;
	}
}
*///commented by Harsha end

template<gpu_size_t BLOCK_SIZE, typename Complex, typename Real>
__global__
void phase_shift_kernel(const Complex* d_in,
                        Complex*       d_out,
                        const Real*    frac,
                        const Real*    delay_rate,
                        gpu_size_t     count,
                        const Real*    fringe,
                        const Real*    dd_fringe,
                        gpu_size_t    fft_size,
                        size_t        fft_cnt,
                        size_t        accum_size,
                        float    update_time)
{
  gpu_size_t thread_idx = threadIdx.x + blockIdx.x*BLOCK_SIZE;
  gpu_size_t thread_count = BLOCK_SIZE * gridDim.x;

    // Iterate over values in the input
  for( gpu_size_t idx=thread_idx; idx<count; idx+=thread_count ) {
    Complex a = d_in[idx];
    const Real pi = 3.141592f;
    float coef;
    size_t ant;
    size_t seq = 0;
    Real coeff_mod;
    ant = (int) idx / (count / 2); 
    seq = (int) idx / fft_size;
    coef = frac[ant] * pi / fft_size;
    if(seq >= accum_size) 
   {
     seq = seq - accum_size;
    }   
  coeff_mod =((idx%fft_size)*(coef - (delay_rate[ant]*pi*update_time*(seq+fft_cnt)/180.f))) + ((fringe[ant]*pi/180.f) - (dd_fringe[ant]*update_time*(pi/180.f)*(seq+fft_cnt))); 
    Complex phase;
    // TODO: Check if pre-computing these and loading from mem is faster
    __sincosf(coeff_mod, &phase.y, &phase.x);

    Complex result;
    result.x = a.x*phase.x + a.y*phase.y;
    result.y = a.y*phase.x - a.x*phase.y;
    d_out[idx] = result;
  }
}



/*//commented by Harsha start
template<typename Complex, typename Real>
void phase_shift_device(const Complex* d_in,
                        Complex*       d_out,
                        gpu_size_t     count,
                        Real           delay,
                        cudaStream_t   stream = 0)
{
	enum { BLOCK_SIZE = 256 };
	dim3 block(BLOCK_SIZE);
	dim3 grid(math::min(math::div_round_up(count,
	                                       (gpu_size_t)BLOCK_SIZE),
	                    (gpu_size_t)CUDA_MAX_GRID_SIZE));
	
	const float two_pi = 2.f*3.141592f;
*///commented by Harsha end
	/*Complex coef;
	coef.x =  cos(two_pi * delay) / count;
	coef.y = -sin(two_pi * delay) / count;
	*/
/*//commented by Harsha start	
  Real coef = two_pi * delay / count;
	phase_shift_kernel<BLOCK_SIZE><<<grid, block, 0, stream>>>(d_in,
	                                                           d_out,
	                                                           coef,
	                                                           count);
#ifdef _DEBUG
	cudaStreamSynchronize(stream);
	cudaError_t error = cudaGetLastError();
	if( error != cudaSuccess ) {
		throw std::runtime_error(std::string("In phase_shift_kernel: ") +
		                         cudaGetErrorString(error));
	}
#endif
}
*///comented by harsha end


template<typename Complex, typename Real>
void phase_shift_device(const Complex* d_in,
                        Complex*       d_out,
                        gpu_size_t     count,
                        const Real*    frac,
                        const Real*    delay_rate,                       
                        const Real*    fringe,
                        const Real*    dd_fringe,
                        float    update_time,
                        gpu_size_t    fft_size,
                        size_t        fft_cnt,
                        size_t        accum_size,
                        cudaStream_t   stream = 0)
{
  enum { BLOCK_SIZE = 128  };
  dim3 block(BLOCK_SIZE);
  dim3 grid(math::min(math::div_round_up(count,
                                         (gpu_size_t)BLOCK_SIZE),
                      (gpu_size_t)CUDA_MAX_GRID_SIZE));

//  const float pi = 3.141592f;
  /*Complex coef;
  coef.x =  cos(two_pi * delay) / count;
  coef.y = -sin(two_pi * delay) / count;
  */
 // Real coef = frac;
 fprintf(stdout, "values of delay at gpu are  %lf %lf %15.12lf %15.12lf\n",frac[0],frac[1],delay_rate[0],delay_rate[1]);
 // Real coef = pi * frac / fft_size;
  phase_shift_kernel<BLOCK_SIZE><<<grid, block, 0, stream>>>(d_in,
                                                             d_out,
                                                             frac,
                                                             delay_rate,
                                                             count,
                                                             fringe,
                                                             dd_fringe,
                                                             fft_size,
                                                             fft_cnt,
                                                             accum_size,
                                                             update_time);
#ifdef _DEBUG
  cudaStreamSynchronize(stream);
  cudaError_t error = cudaGetLastError();
  if( error != cudaSuccess ) {
    throw std::runtime_error(std::string("In phase_shift_kernel: ") +
                             cudaGetErrorString(error));
  }
#endif
}

